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^ ; ABSTRACT 



We determine an empirical dense matter equation of state from a heterogeneous dataset 
of six neutron stars: three type I X-ray bursters with photospheric radius expansion, studied by 
Ozel et al., and three transient low-mass X-ray binaries. We critically assess the mass and radius 
O ' determinations from the X-ray burst sources and show explicitly how systematic uncertainties, 

o 
o 



lO ' such as the photospheric radius at touchdown, affect the most probable masses and radii. We 

^^ ■ introduce a parameterized equation of state and use a Markov Chain Monte Cai^lo algorithm 

within a Bayesian framework to determine nuclear parameters such as the incompressibility 

and the density dependence of the bulk symmetry energy. Using this framework we show, for 

k> ', the first time, that these parameters, predicted solely on the basis of astrophysical observations, 

H ' alllie in ranges expected from nuclear systematics and laboratory experiments. We find signifi- 

ed ' 

cant constraints on the mass-radius relation for neutron stars, and hence on the pressure-density 

relation of dense matter. The predicted symmetry energy and the equation of state near the 
saturation density are soft, resulting in relatively small neutron star radii around 11-12 km for 
M = 1.4 Mq. The predicted equation of state stiffens at higher densities, however, and our pre- 
ferred model for X-ray bursts suggests that the neutron star maximum mass is relatively large, 
1.9-2.2 Mq. Our results imply that several commonly used equations of state are inconsistent 
with observations. 

Subject headings: dense matter — stars: neutron — X-rays: binaries — X-rays: bursts 



1. Introduction 



The masses and radii of neutron stars are determined by the pressure-energy density relation ( equation 
of sta t e; EOS hereafter) of cold dense matter using the familiar Tolman-Oppenheimer-Volkov (TOV dTohnan 
19391 : lOppenheimer & VoUcofil 1 19391) relativistic stellar structure equations. Within tens of seconds after 
birth, the neutron star is cold (meaning temperature much less than Fermi energy) and deleptonized (meaning 
in beta equilibrium with no trapped neutrinos); as a result, for a given EOS the mass and radius of the star 
depend only on the central density. Inversion of the structure equations, given simultaneous mass and radius 
measurements, can therefore constrain the pressure-density relation, although the quality of the constraints 
are very sensitive to the mass and radius measurement uncertainties. 

A host of observable phenomena and experimental information is becoming available (for a recent 
review, see 



Lattimer & PrakashI |2007|) : on the observational side, pulsar timing, thermal emission from 



cooling neutron stars, surface explosions, and gravity wave emissions are some promising areas for mea- 
surements of mass and radius; on the experimental side, heavy-ion collisions, giant dipole resonances, and 
parity-violating electron scattering are some promising techniques for measuring the density dependence of 
the pressure of nuclear matter. These efforts ai^e complementary, and determining an EOS from the many 
disparate measurements is a challenging task. 

This paper consists of two parts. In the first, we consider simultaneous mass and radius information 
from astrophysical observations of X-ray bursts and thermal emissions from quiescent low-mass X-ray bi- 
naries. We critically asses s the mass and radius constraints determined from X-ray bursts that may reac h 
the Eddington limit (§ E]) (|vanParadijslll979l . Il982l : IPaczynski & AndersoJ 19861 : Ivan Paradijs et al.lll99d) . 
Heretofore, these bursts have been interpreted assuming that the photospheri c radius is equal to the stellar 
radius at "touchdown", when the effective temperatur e reaches a maximum (|Ozell 120061 : lOzel et al.ll2009l : 
Gilver et al.ll2010l : lGuver et al.ll2010bl : IOzel et al.ll2010l) . If true, rather stringent constraints on mass and ra- 
dius are predicted. Using the most probable vali ies for the observe d flux and angular emission area, however, 
does not result in real- valued masses and radii. lOzel et al.l (120101) argues for rejecting much of the observed 
phase space on this basis and thus obtains tight constraints on masses and radii. We show that this model 
is not internally consistent, but that consistency can be regained by relaxing the assumption concerning the 
effective photospheric radius. With this modification, the inferred neutron star radii moderately increase 
and con fidence intervals for predicted masses and rad i i become substantially larger than th ose previously 
quoted Jozellbood : bzel et al.ll2()()9l : lOuver et al]l201ol : lGuver et al.ll2010bl : IOzel et allboiol) . 



We then discuss (§ O mass and radius estimates from thermal emission from quiescent low-mass X- 
ray binaries (LMXBs hereafter). Neutron stars in quiescent low-mass X-ray binaries may be used to obtain 
angular emission areas, and thereby mass and radius information, from their quiescent thermal emission. In 
total, we have mass and radius constraints for six neutron stars: three bursting sources with photospheric 
radius expansion bursts and three quiescent neutron stai^ transients. 

While not all of the uncertainties involved in constraining the masses and radii of neutron stars are under 
control, it is important to quantify the constraints on the EOS which are implied by the observations. In the 
second part of this paper (§ |4l), we use a Bayesian analysis to combine these mass and radius constraints 



to determine empirical pressure-density and neutron star mass-radius relations. We include constraints to 
the EOS from causality, the observed minimum value of the neutron star maximum mass, and the observed 
maximum pulsar spin frequency. We employ a parameterized EOS (§ 14.21 ) that is compatible with laboratory 
measurements, and use it with a Markov Chain Monte Carlo (MC) algorithm to jointly fit these mass-radius 
constraints (§ 14.31 ). We show, for the first time, using astrophysical observations alone, that values of the 
nuclear parameters such as the incompressibility K, bulk symmetry energy 5,,, and the density dependence 
of the symmetry energy y, all lie in ranges expected from nuclear systematics and laboratory experiments. 
Furthermore, our results imply that the most likely neutron star radius is relatively small, of order 1 1- 
12 km for neutron stars with masses near 1.4 M©, so that the predicted EOS is relatively soft in the density 
range 1-3 times the nuclear saturation density. Implications for the maximum mass sensitively depend on 
assumptions concerning X-ray burst models and ai^e discussed in § 14.41 In §5 we discuss astrophysical and 
nuclear physical implications of our results, and compare our methods and results with other studies. 



2. Mass and radius constraints from photospheric radius expansion bursts 



Although a few dozen neutron star masses have been determi ned very accurately (to w ithin a few 
percent) in binaries containing pulsai^s (for a recent compilation, see iLattimer & PrakashI 120051) . no radius 
information is available for these systems. For systems in which the neutron star accretes matter from 
a nearby companion, nuclear processes in the crust and envelope of the neutron star provide additional 
obseryables that can be used to constrain its mass and radius. Our discussion here initially follows that of 
Qzell (120061) . and is provided to establish a framework for the subsequent assessment of this method in §§ 12.11 
and[ 



Type I X-ray bursts are the result of thermally unstable helium (or in some cases, hyd rogen) ignition 
in the accreted envelope of a neutron star (for a review, see IStrohmayer & Bildstenll2004l) . The ignition 
generates a thermonuclear explosion that is observed as an X-ray burst with a rapid rise (~ 1 s) followed by a 
slower decay (~ 10-100 s). With the discovery of significant spectral softening during some bursts it quickly 
became apparent that significant radial expansion of the photosphere can occur during powerful X-ray bursts: 
if the burst is sufficiently luminous, radiation pressure drives the photosphere outwards to larger radii, in 
some cases subs tantially so. About 20 % of X-ray bursts show evidence for photospheric radius expansion 
(hereafter PRE; iGallowav et all 120081) . The inference that radius expansion occurred spurred many (e.g.. 



Paczvriski 1983 



Ebisuzaki et al. 



19831) to construct models of extended, radiation-dominated neutron star 



envelopes. There is widespread agreement from these calculations that during a radius expansion burst, 
the flux at the photosphere approaches (to within a few percent ) the Eddington value. The convective zone 
is expected to reach the photosphere during a powerful burst dWoosley et al.ll2004l : [Weinberg et al.ll2005l) 
thereby polluting the accreted material with heavier nuclei synthesized during the burst. 

During a typical PRE burst, the flux rapidly increases, peaks, and then decreases. While the flux f oo 
is near maximum, the blackbody temperature Tbb.oo at first decreases, then increases to a maximum before 
decreasing again. (The co subscript indicates that the quantities are observed at the Earth and therefore 



differ from their value in a reference frame local to the emission.) During this time, the normalized area 
^oo/7"bbco increases to a maximum value and then decreases as rbb,oo increases. The point at which Tbb.oo 
reaches a maximum (and the normalized area typically stops decreasing and becomes constant) is thought to 
be when the photosphere "touches down" at the stellar radius. We shall denote the observed flux, measured 
at this time, as f td.oo- Following this point, the angular emission area remains roughly constant as both 
Fco and Tbb.oo slowly decrease. At the peak of the expansion, the low blackbody temperature puts much 
of the emitted flux below the ban dpass for many X-ray instruments; in extreme cas e s this can create the 
appearance of a "precursor" burst (IHoffman et al.lll978uin 't Zand & Weinbergll20ld) . iDamen et al.l (Il990l) 
argued that to minimize systematic errors in determining the Eddington flux, the measurement should be 
made at touchdown, when the the temperature is at a maximum and the photosphere has presumably just 
retreated to near the quiescent stellar radius. 

At touchdown, the observed flux is expected to equal the local, i.e., as measured in a frame co-moving 
with the photosphere, Eddington value 



TD,oo 



GMc 



^1 - 2/3(rph) = f Edd ^1 - 2y6(rph). 



(1) 



HereyS(r} = GMI{rc^), k is the opacity, and rph is the photospheric radius at the time this flux is evaluated. 
Ozell (120061) argued that at touchdown, i.e., when Tbb reaches a maximum, rph - R. For clarity, when we 
refer to the Eddington flux in the remainder of the paper, we shall mean Fgdd = GMc/{kD^). This definition 
is independent of the stellar radius, and is a true limiting flux: Ftd,oo < ^Edd, with equality holding if 
rph » R. Finally, for Thomson scattering in a hydrogen-helium plasma, k « 0.2(1 -i- X) cm^ g~\ where X is 
the mass fraction of hydrogen. 



Galloway et al.l (120081) examined a large sample of PRE bursts. They found that in all cases for which 



the inclination was not extremely high, the touchdown flux was within a factor of 1.6 of the peak flux, 
consistent with it occumn g when the photosphere had retreated. The source EXO 748-676, which wa s 
used in previous analysis (iQzell l2006l) because of its claimed redshift measurement (ICottam et al.l 120021) . 
has a high inclination (it is a "dipper"); as a result, the observed touchdown flux may be obscured by the 
disk, and indeed f td.oo is much less than the observed maximum flux for this source. The distance to 
EXO 0748-676 is also not well known, and for these reasons we omit this source in our analysis below. 

In the latter part of the burst the ratio FoolT^-^ ^ is observed to be roughly constant. This allows one to 
define a normaUzed angular surface area. 



bb,oo 



-in^. 



(1 - 2y6) 



-1 



(2) 



Here Fco and Tbb.oo are the flux and blackbody temperature as measured by a distant observer in the late 
phase of the burst, when A is roughly constant, /? = GM/{Rc^) is the compactness, D is the distance to the 
source, and fc = T],b/Tss is a color correction factor that accounts for the departure of the spectrum from a 
blackbody. 



If a distance to the star can be estimated with sufficient precision, the observed touchdown flux (when 
the blackbody temperature reaches a maximum following the PRE) and the inferred apparent angular area 
measured late in the burst can be converted into an estimate of the stellar mass and radius. It has been 



Ozel et al. 



Ozel et al. 



clairned that 1-cr uncer tainties of less than 10% in both mass and radius are pos sible ( Qzelll2006l : 
20091: lOzel et al.ll2010l) . Systematic uncertainties not included in the analyses of lOzell (120061) and 
(120091 ) include the composition of the accreted material and th e effects of the neutron star atmosp here on 
the spectral shape, although these are included in later analyses (IGiiver et alJl2010l : IOzel et alJl2010l ). A key 
question is whether the value of yS(rph) in equation ([T|l, which is determined by rph, is the same as the value 
of yS in equation (0. In the following discussion, we shall first assume that these y6 values are the same, 
i.e., Tph = R, as would be the case if the photosphere has indeed just retreated to the stellai^ radius. We will 
demonstrate that for EXO 1745-248, 4U 1608-522, and 4U 1820-30 the most probable observed values of 
A and FEdd,oo do not lead to real-valued solutions for M and R. Indeed. lOzel et all (120101) notes that forcing a 
real-valued solution for mass and radius by Monte Carlo sampling within the probability distributions of the 
observables Ftd,oo, A, and D reduces the uncertainties in the mass and radius to values smaller than those of 
the measurements. We shall explore in § 12.21 a diff"erent interpretation, which is that r^^ > R when Fxd,oo is 
evaluated, so that f td.oo ^ f Edd- We show that relaxing the assumption rph = R generally yields real-valued 
solutions for M and R for the most probable values of the observables. 

We combine the observed quantities Ftd.co and A, along with a measurement of the distance D and 
theoretical estimates of fc and k, into two parameters, 

Ftd,oo kD 



a 



7 = 



FtD,ooK 



(3) 



(4) 



If we then make the assumption, following lOzell (120060 . that f td,oo = ^Edd V^"^^^ ^^ equation ([T]), we find 



a 



Al - 2yS), 

R 



y6(l-2;S)3/2' 

Solving eq. <^ for jS, we then solve eq. ^ for the radius and mass: 



R 

M 



1 1 r- 

4 4 



8a, 



arVi-2jS, 

G 



(5) 
(6) 



(7) 
(8) 
(9) 



Note that y is independent of D. An important consequence of eq. (O is that for both M and R to be real, 
we must have a < 1/8. Since a is determined, however, from the observables Fjd.oo, A, and D, as well 
as the estimated parameters k and fc, the inferred value of a does not necessarily satisfy this mathematical 
limit. This condition serves as a check on the validity of the assumptions made in the modeling of the radius 
expansion bursts. 



2.1. Monte Carlo analysis of photospheric radius expansion burst data 



Observational inform ation for three Type-I X-ray bursters with PRE bursts (iQzel et al.ll2009l : lGuver et al 



2010l : lGuver et al.ll2010br) is presented in Table [T] along with associated uncertainties. We have reexamined 
the uncertainties for each object as described in Appendix |Al in some cases our values for the uncertainties 
differ from those used previously. 

We first fix fc ^ 1.4 and X ^ as was done in previous work for EXO 1745-248 Jozel et allboooh . We 
observe in this case from Table [1] that no real- valued solution of equations Q-® is possible for any of the 
sources because a > 1/8 for the central values of the observables. Moreover, it is not consistent with X-ray 
burst models to fix the color-correction factor to /^ = 1.4 and the hydrogen abundance to X - 0, and indeed 
these restrictions w ere relaxed in subseque nt analyses of 4U 1 608-522 and 4U 1820-30 (iGiiver et al.ll2010l : 



Gilver et al.ll2010bl) . Atmosphere m odels (IMadei et al. 



consistent with earUer calculations (ILondon et al.lll98q : lEbisuzaki & Nakamuralll988h . The largest values 



20041) suggest a range 1 33 < fr < 1.81; these are 



of fc are approached as the flux approaches the Eddington limit. In the tail of the burst, the flux is lower and 
fc does not vary strongly. We choose to select the value of fc in our Monte Carlo analysis from a boxcar 
distribution centered at /^ = 1 .4 with an uncertai nty of 0.07. This selection is compa rable to the boxcar 
distribution with/c = 1.35 ±0.05 used previously (iGiiver et al.ll2010l : lGuver et al.ll2010br) . (In the following, 
we denote a boxcar uncertainty inZ as AX and a Gaussian uncertainty with crx-) In addition, as described in 
Appendix lA. 1 l and lA.2[ we find insufficient information from the observations to exclude any possible values 
of X for EXO 17 45-248 and 4U 1608-522, so we choose a uniform distribution with < X < 0.7, which is 
also assumed by IGiiver et al.l fcOlOi ): IGiiver et al.l (l2010bh . The source 4U 1820-30 is an ultracompact and 
accretes He-rich fuel; see Appendix IA.3I 

If values of fc and X or the observables Fjd.co, A, or D are selected at random from their respective 
probability distributions, some combinations will satisfy the condition a < 1/8. In this way, using only these 
acceptable combinations, most probable values a and y can be established using Monte Cai^lo sampling. 
(Here and below, we define X to be the most probable value of X obtained from a MC simulation, and 
the uncertainties on this quantity ai^e determined by selecting the region surrounding the most probable 
value that includes 68% of the total MC weight.) Values for M and R, and their uncertainties, can then be 
determined from a and y using equations (O-®. But the fraction of accepted realizations is then very small 
as shown in Table |2] The entries in the first group of this table give the most probable values a, y, and ^oo 
using Tph = R and values for Ftd,oo,A, and D from their probability distributions summarized in Table [T] 
and, for fc and X, from previous discussion. For each quantity, the uncertainty is determined by creating a 
histogram for the indicated MC realizations, sorting the bins by decreasing weight, and selecting the range 



which encloses 68% of the sum of all bins. We identify 7?oo = [FooK^tto-T^^^)] ' = /?/ y 1 - 2/3 = ay 
and impose no a priori restriction for the value of a. The most probable value of a, a, is observed to be 
approximately a factor 1 -i- X larger than the value of a in Table \T\ where X is the average value of X in 
its assumed range, i.e. 1.35 for EXO 1745-248 and 4U 1608-522, and 1.0 for 4U 1820-30. For all three 
sources, a > 1/8 by several standard deviations. 



Table 1. Observational values for Type I X-ray burst sources used in this paper; see Appendix lAl for 
details about how the uncertainties were determined. 



Quantity 



EXO 1745-248 4U 1608-522 



4U 1820-30 



D (kpc) 

A (km^ kpc"2) 

^TD,oo (10"^ ergs cm~^ s~^) 

y (kmf 

Roo = ay (km)'^ 



6.3 ± 0.6 
1.17 ±0.13 

6.25 ± 0.2 

0.131 ±0.017 

101.7 ±11.8 

13.4 ±2.0 



5.8 ± 2.0^ 
3.246 ± 0.024 0.9198 

15.41 ±0.65 
0.179 ±0.062 

114.5 ±4.9 

20.5 ±7.1 



8.2 ± 0.7'' 

0.0186 

5.39 ±0.12 

0.166 ±0.015 

92.7 ± 2.8 

15.4 ±1.4 



References. 



dzelll2006l : lOzel et al.ll2009l : iGuver et al.ll20inl : iGuver et al.ll2010b 



^The distance distribution for 4U 1608-522 was cut off below 3.9 kpc. 



''The distance uncertainty for 4U 1820-30 was approximated by IGuver et al. 
(|2010bh as a boxcar with halfwidth 1.4 kpc. 



'^Assuming fc- 1 .4 and X = with A/^ = AX - 0. Errors are computed assuming 
uncorrected Gaussian errors for D, Fyd oo and A. 



Table 2. Comparison of Monte Carlo realizations for PRE bursts employ ing data from Tab l e [D In the 

thi rd gro u p, we use the values and distributions of Fjd, A, D, fc, and X from lOzel et al.1 (l2009h : lGuver et al 

(120101 ): iGiiver et al.1 (l2010bh . For all other entries, we use < X < 0.7 for the hydrogen mass fraction, 

except for 4U 1820-30, for which X -0 and we use 1.33 < fc < 1.41 for the color correction factor. 



EXO 1745-248 4U 1608-522 



4U 1820-30 



TD.oo 



^Edd[l - 2yS(/?)]'/2. Q, um-estricted 



a 



y (km) 
.^oo (km) 



165+°°'^^ 
70 4+''^-2 

'"•^-14.7 

13 2+'-^ 



A 99Q+0.054 
"•^^^-0.038 

20.5!^:^ 



80-91.^56 
15.3!}- 



i^TD,oo - /^Edd[l - 2piR)]^'^; or < 1/8 restriction 



a 


A 199+0.003 
"■^^^-0.007 


123+0003 

"■ ^^-'-0.005 


Q 123+0.003 

"•^^■'-0.004 


f(km) 


107 9+''*-2 

^"'•^-14.2 


127 1+^° 


109.3+5-] 


/?co(km) 


19 7+1-6 


1 c 9+0.85 
^■^■^-0.81 


13.2!?;f 


Points accepted 


4.4% 


0.24 % 


0.44% 



Ftd,oo = fEdd[l - 2p{R)]^'^;a < 1/8; Original inputs from Ozel et al. 



a 

y(km) 
^oo(km) 


A 122+0003 

"•^^ -0.003 

109.0!^-2 
1 3 3+0.4 


A 199+0.003 
"■ ^^^-0.003 

113.1!^;° 
13.8!°:^ 


124+0001 

"•^^ -0.001 

104.0!}-^ 

1 9 Q+0.2 
^^•^-0.2 


Points accepted 


13% 


0.18% 


1.5x10"^% 



TD.oo 



FEdd; a < 3~^'^ = 0.192 restriction 



a 

7 (km) 

^oo(km) 


A 1 f:5+0.026 

79 0+i''-2 

'^■"-12.4 

13 1+'-^ 


97.0!}^-4 
15.3!2j 


A 1^:7+0.016 

"■^'^'-0.016 

^^■"-10.6 
15.0!};^ 


Points accepted 


65% 


15% 


91 % 



FEddU - 2piR)]^'^ < Ftd.co < fEdd;^^ +p^ <0 restriction 



a 

7 (km) 

^00 (km) 

Points accepted 


0-143!°:°?^ 
87.1_+}^i 
13.0!};^ 
32% 


155+0018 

103.0+}^-^ 
15 3+2-3 ■ 

5.8% 


0.158!°:°}^ 

104.8!^,^o 
15.0!};^ 

37% 



In the second group of Table [2l we repeat the analysis, but this time only accept realizations for which 
a < 1/8. Note that at most 4% of the realizations are accepted. This implies, at a minimum, that the assump- 
tions in this model are incomplete. The most probable value for a decreases and now lies approximately 
1-cr below the value 1/8, and the size of confidence interval for a is now significantly reduced compared 
to the unrestricted case, by factors of 4-18. The net effect is to preselect the value 1/8 for a irrespective 
of the observed values of f td.oo, A, and D. Another consequence of the selective rejection of most of the 
Monte Cai^lo realizations is to increase the most probable value of y and decreases the r nost probab l e valu e 
of /?oo, and to greatly decrease their uncertainties; the latter result was already noted by lOzel et all (120101) . 
In addition, the r e sulting values of y and Rqq for the diff erent sources are also "herded" into similar values 
(lOzel et al-lboogl : lOuver et allboiol : lOuver et allboiObh . The third group summarizes the MC results ob- 
tained using the exact same observational inputs found in those references. Note that for 4U 1820-30, only 
about 1 realization out of 100 million is now accepted! The larger acceptance rate for EXO 1745-248 is 
due to the assumption that X - which gives the smallest possible a. Analyses of X-ray bursts with this 
model result in certain values for M and R that are neai^ly independent of a and y, and are therefore almost 
completely independent of the observables Fjdoo,^, and D, which seems unrealistic. The fact that only a 
small fraction of realizations are accepted may indicate that the underlying model contains systematic errors 
or is unphysical, a point we shall explore in the next section. 

Figure [U displays the Monte Carlo probability distributions for the masses and radii of EXO 1745- 
248, 4U 1608-522, and 4U 1820-30 as determined from a and y following equations ([T])-®, and with 
the restriction a < 1/8. There are two pe aks in the distributio ns because of the quadratic in equation ([7]). 
We also impose thatyS < 1/2.94 ~ 0.34 (lGlendennindll992l'). whi ch is based on requiring a subluminal 
sound speed throughout the star (see also iLattimer & PrakashI 120070 . This causes the rejection of a small 
fraction of realizations from the higher- redshift region. For all three sources, our confidence intervals are 
larger than computed in previous works (|Qzel et al.ll2009l : lGuver et al.ll2010l : lGuver et alEOlObT) . The origin 
of this distinction is because we use larger variations fc and, for EXO 1745-248, in X, and we employ a 
different uncertainty in D for 4U 1820-30. A detailed discussion of the treatment of the uncertainties in the 
observables for each object is given in the Appendix. 

When using a Monte Carl o scheme to determine uncertainties, \y e find that the two err or ra nges in Fig. 



[T] are e qually populated, while lOzel et all (|2009h . lOuver et all (|2010h . lOuver et all pOlObI '). and lOzel et al. 



20101 '). transforming probabilities using a Jacobian, indicate a much higher probability for the higher redshift 
solution. This difference is caused by a numerical error in their Jacobian (Ozel, priv. comm.); when the 
correct Jacobian is used, both integration methods agree. 



2.2. Alternate interpretations of the location of the photospheric radius at touchdown 



The lack of real-valued solutions for M and R for the most probable values of the observables, and 
indeed, the rejection of the vast majority of Monte Carlo realizations, motivates us to consider another 
possibihty; namely, that at "touchdown" the photosphere is still extended. To explore this situation, we first 
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Fig. 1. — Mass-radius probability distributions for Type I X-ray bursts assuming that the photospheric radius 
and the stellar radius are identical. The causal limit yS = 1/2.94 is indicated with a dashed line. These plots 
con^espond to the results shown in the second group in Table[2l The solid curves indicate the 68% and 95% 
confidence boundaries while the shading level reflects the relative probabilities. All distributions, P,, are 
normalized so that J f , dM dR = I. 



11 



consider the extreme possibility that rph » 7? at the point identified as touchdown from the maximum in 
Tbb.oo- In this case ^ td,oo = f Edd and is independent of the stellar radius R. With this assumption, a (eq. lO) 
and y (eq. JU) are now related to p and R via 

a - /3yJ\-2p, (10) 

R 

7 = r,,. ... • (11) 



Ai - m 



Defining a quantity Q - cos (1 - 54q' ), the expressions for the compactness, radius, and mass are then 



R 

M 



- [l + V3 sin(0/3) - cos(6i/3)] , 
^[l + 2cos(6i/3)]. 



ay 



4r^ 



1,2- 



— ot y, 
G ^ 



(12) 

(13) 
(14) 

(15) 



Obviously, M is real for all a, y > 0. For a < 3^^''^ ^ 0.192, Q is real and there are 3 real roots foryS and R. 
Only two of these, < ySi < 1/3 and 1/3 < yS2 < 1/2, are physically meaningful: the other root is negative. 
When a > 3"^^^ there is one real root for j6 and two imaginary ones, but the real root is negative. 

From the top group in Table El we can see that both 4U 1820-30 and EXO 1745-248 have most- 
probable values of a that satisfy the constraint for positive real values of yS and R, and 4U 1608-522 lies less 
than l-tr above this limit. Therefore, we now find that a much larger fraction of the Monte Carlo realizations 
are accepted when selecting values for Ftd,oo, A, D, f^, and X from their probability distributions. This is 
shown in the fourth group of Table |2l for which we use the model given in equations (ITOll-dTTI) and impose 
the restriction a < 0.192. In contrast to the case in which rph = R, the uncertainties in a and y are not as 
strongly diminished. Moreover, the values for a are no longer nearly the same for the three sources. While 
in principle each accepted realization results in two M, R values, one corresponding to ySi and the other to 
P2, nearly all theyS2 realizations are rejected on the basis of causality when/?2 > ^ 12.9 A. Figured displays 
the probability distributions for the fourth group of runs listed in Table |2l as well as their 68% and 95% 
contours. T he eqor contours are larger than those shown in Fig. dH) and conside rably larger than determined 
previously dozellbood : bzel et al.lEo09l: lOuver et aPboidbuver et allboiObl) . 



It is perhaps unphysical to make the extreme assumption rph » R when T^s reaches a maximum. But 
because of the relative consistency of the solutions and the much larger fraction of MC points accepted in 
this case, it seems reasonable to consider the possibility that at touchdown, that is, when Tbb.oo reaches a 
maximum and the normalization reaches a minimum, the photospheric radius may in fact be larger than the 
quiescent stellar radius. We can write a = yS -^^T^^^ ^/1~^^^, where ySph = GM/(rphC^) < /3. Defining the 
quantity h = 27?/rph, the quantities a and y now become 

f Edd kD 



a 



3f2 



Ac^f- 



■./3yJ\-2/3yJl-hf3, 



R 



FEddK p{l-2p)^Jl-h/3 



(16) 



(17) 
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Fig. 2. — Mass-radius probability distributions for Type I X-ray bursts assuming that rph » R. The causal 
limit yS = 1/2.94 is indicated with a dashed line. These plots correspond to the results shown in the fourth 
group in Table|2l The solid curves indicate the 68% and 95% confidence boundaries while the shading level 
reflects the relative probabiUties. All distributions, P,, are normalized so that j Pj dM dR - I. 
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Note that /?„ = ay remains unchanged. Equation (1161) is a quartic equation for/3. Two positive real solutions 
for y6 exist when q^ + p^ < 0, where 



P = ^1'^' 



1 

lAh] ' 



-Ti^dii^^^""'!' 



^ ^ ^'^ 



3(2 + /?)^ 
16/j 



(18) 
(19) 
(20) 



otherwise the only real solution is negative. Defining the quantities 

,3/2" 



= COS 



-1 



V 

w 



b 



= 2^f^cos(e/3), 

^ V2v - 2a , 

-\2v + 4a± — 
w 



(2 + hf 
8/i2 



1- 



2 + h 



8/7 



(21) 

(22) 
(23) 

(24) 
(25) 



the two solutions for /3 with values in the range < y6 < 0.5 are 

2 + h w - Vz+ 

where the sign has the same sense in the two occurrences. 

Figure [3] displays the resulting probability distributions assuming that h is uniformly distributed in the 
range < h < 2, i.e., R < r^h < oo. The eiTor contours are much larger than those determined in the h = 2 
(rph = R) case and slightly larger than those in the h = (rph » R) cases. Most of the solutions for f3+ 
are rejected because they lead to acausal combinations of M and R. The fifth group in Table |2] presents the 
associated values of a, y, and /?£». These results are nearly identical to those of Figure |2] because small 
values of the photospheric radius are strongly disfavored by the requirement that the masses and radii are 
real-valued. We will therefore use the rph » R probability distributions in M and R for the Bayesian 
estimation of the EOS in Section IH and note that these results will be essentially identical to what one 
would obtain without assuming a particular value for the radius of the photosphere, as long as it is not a 
priori assumed to be equal to R. 

Using the number of accepted MC realizations as a guide, we can also obtain an estimate for the lower 
limit of the location of the photosphere at touchdown. As is common in Bayesian analysis, we use the ratio 
0.1 as a guide; models for which the fraction of accepted realizations is less than 0.1 are rejected. This 
imphes that rph > 5.07? for 4U 1608-522, rph > 1.1/? for EXO 1745-248, and rph > 1.4/? for 4U 1820-30. 
This analysis appears to strongly disfavor an interpretation of the touchdown radius rph = R and we find 
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it likely that the photospheric radius is extended at touchdown. At the same time, there is little difference 
between the results assuming /j = or a uniform range < h < 2. It would then seem justified that we, in 
our remaining analysis, reject the interpretation rph = R, and make use of the rph » 7?, /j = interpretation. 
Nonetheless, we will also include results for the rph = R interpretation in the work below for comparison 
to previous studies. We will show that some aspects of the EOS constraints are dependent on whether one 
assumes that rph = R or rph » R. 



3. Mass and radius constraints from thermal spectra 



The masses and radii of isolated neutron stars or ones in transient LMXBs can be inferred from spectral 
modeling if their distances are accurately determined. Relatively accurate distances are known for several 
accreting neutron star transients located in globular clusters. Although the uncertainties for any individual 
source are large, it is productive to use the ensemble of observations to improve constraints on the dense mat- 
ter equation of state. In addition, future mass and radius measurements from thermal sources are expected 
to tighten such constraints further. 

Many neutron stars are in transients, for which the accretion of matter proceeds intermittently, with 
episodes of accretion separated by long periods of q uiescence. While the ne utron star accretes, compres- 
sion of matter in the crust induces nuclear reactions (IHaensel & Zduniklll990l) that rel ease heat. When th e 



accretion ceases, the heated crust cools, resulting in an observable thermal luminosit y (IBrown et al.l 



Because the timescale for heavier nuclei to sink below the photosphere is short (~ 10 s: lBildsten et al. 



1998 ^. 



1992b, 



and these systems show no evidence, such as pulsations or cyclotron spectral features, for a significant 
magnetic field, the spectra can be fitted with well-understoo d unmagnetized hydrogen atmosphere spectra 
(|Zavlin et al.lll996l : iRutledge et al.lll999l : iHeinke et al.ll2006r) . As a result, the observed X-ra y spectra can 



be us ed to reliably infer an apparent angular emitting area, and, possibly, the surface gravity (IHeinke et al. 



20061) . Such objects with accurately determined distances (such as those in globular clusters) can be used to 
estimate masses and radii. 

The spec tra of isolated cooli ng neutron stars with well-determined distances, such as RX Jl 856-3754 
discovered by IWalter et al.l (119961 ). have a much larger signal-to-noise than those of neutron star transients. 
The interpretation of their spectra is complicated, ho wever, by the poten t ially s trong magnetic field. The 
distance to RX J 1856-3754 has been controversial. IWalter & Lattimed (120021) gave a pai'al lax distance 
D = 117 + 12 pc based on HST observations from 1996-98, b ut Ivan Kerkwijk & KaplanI (120071) later found 



Walter et al. 



2010l) reanalyzed the new data 



D = 161 ± 16 pc using new HST data from 2002-4. Recently, 

and found D - 122 ± 13 pc, in good agreement with their older estimate. A multi-wavelength spectral fit 



by iPons et al.l (120021) for non-magnetic heavy element atmospheres (which they argued gave the best fits to 
the X-ray and optica l spectra of RX J1845-3754) yielded Ro^/D ^ 0.13 ± 0.01 km pc^^ and 0.3 < z < 0.4. 
Burwitz et all (|200ll ) argued, however, that these models predict spectral f eatures such as absorption edges 
that are not apparent in the data. In addition, the ex istence of a bow shock (Ivan Kerkwijk & Kulkamill200ll ) 



and the detection of pulsations in the X-ray flux (JTiengo & Mereghettil l2007l ) and their period derivative 
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Fig. 3. — Mass-radius probability distributions for Type I X-ray bursts assuming a uniform distribution in 
h = 2/?/rph. The shadings and lines have the same meaning as in Fig.|2l 
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(Ivan Kerkwijk & KaplanI 120081) imply this star has a magnetic field of order lO'^G. Nevertheless, there 
could be many reasons why spectral features are washed-out from a rotating, highly magnetized, object. 
A two-temperature blackbody model, in which the X-ray flux is primarily emitted by a high-temperature 
region and the optical flux by a low-t emperature region, gives Rco/D - 0.14kmpc~\ approximately the 
same value d etermined bv IPons et all ( 2002J) and iBurwitz et al.1 (120031) . Motivated by the large distance 



determi ned bvlvan Kerkwiik & KaplanI (|2007h . which implied radii incompatible with neutron star equations 
of state, IHo et al.l (|2007h developed a magnetized neutron star atmosphere model with a condensed surface 
and a trace amount of H remaining in the atmosphere. This model resulted in Rco/D - 0.12kmpc~^. 
The origin of the trace H in the atmosphere is u nknown, and its m ass must be finely-tuned to explain 
the magnitude o f the optical flux. T he analysis of 



Pons et al. 



(120021) coupled with the confirmed distance 
determination of [Walter et all (1201 Of) yields M = 1 .7 + 0.3 M© and /? = 1 1 .5 ± 1 .2 km, but this model neither 
accounts for a lai^ge magnetic field or the observed lack of spectral lines. Because of this, we do not include 
this source in our baseline results although we have performed a sepai^ate analysis to determine if our results 
would be appreciably affected. 

Some of the str ongest constraints on properties of transient neutron stars ai^e for X7 in the glob- 
ular c luster 47 Tuc (Heinke et al.ll2006l) . and for the neutron stars in the globular clusters o) Cen and 



M13 dWebb & Barret 



20071). The 99% contours which define the M and R values infeiTcd by the atmo- 



Heinke et al. 



sphere models in Fig. 8 of IWebb & Barretl (|2007l) and the 68%, 90%, and 99% contours inferred by the 
atmosphere models in Fig. 2 of 



(I2OO6I) are shown in Fig. |4l For use in constraining the EOS, 
we need to create probability distributions for these observations in the M, R plane. While one can roughly 
represent the published constraints on M and R in terms of a single constraint on Rco, this is accurate only 
for smaller masses. To construct a better representation of the data, we first note that the likelihood contours 
are nearly perpendicular to lines emanating from the origin of the M, R plane. To construct a probability 
distribution, we assume that the probability is distributed as a Gaussian with a width determined by the 
spacing of the contours across a line going through the origin. Figure H] shows a comparison between the 
resulting probability distributions used in this work and contours obtained in the original references. 



4. Application of statistical methods to constraints on the EOS 

As is clear from the preceding discussion, not all of the uncertainties involved in constraining the 
masses and radii of neutron stars are under control. Nevertheless, it is interesting to understand what these 
observations may imply for the EOS. Furthermore, it is important to quantify their implications for the EOS 
in order to motivate future observational work that will reduce these uncertainties. 

In this section, we apply a Bayesian analysis to the data described above. We will first briefly review the 
formalism (§ 14. IK and then develop (§ 14.21 ) a general parameterization of the EOS. The output probability 
distribution for the EOS parameters and the EOS itself are presented in § 14.31 The most probable masses 
and radii are presented in § 14.41 Even though the results of § |2] strongly suggest that rph = /? is disfavored, 
we include results for both rph = R and rph » /? to simplify comparison to previous studies. 
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Fig. 4. — (Upper left panel) The solid curve is the 99% probability contour in the mass-radius plane for the 
LMXRB in M13 from I Webb & Barrel! (120071) . The red shading shows the probability distribution adopted in 
this paper. ( Upper right panel) T he various lines are the 68%, 90%, and 99% probability contours for X7 in 
47 Tuc from lHeinke et al.l (120061 ). The red shading shows the probability distribution adopted in this paper. 
(Lower left panel) The same as the upper left panel except for oj Cen. All distributions, P,, are normalized 
sothat JP;<iM<i/? = 1. 



4.1. Bayesian analysis 



Bayes theorem can be formulated as (see, e.g. jGrinstead & Snellll 19971) 



P{D\M)P{M) 
P{M\D) = '^^^^' , (27) 

where P(M) is the prior probability of the model At without any information from the data D, P{D) is the 
prior probability of the data D, P{D\M) is the conditional probability of the data D given the model M, and 
P{M\D) is the conditional probability of the model At given the data D. This latter quantity, P{M\D) is 
what we want to obtain, namely, the probability that a given model is correct given the data. 

For many non-overlapping models At, which exhaust the total model space At, this relation can be 
rewritten 

^^^'1^^ - i:,p(DiAt,)P(At,) • ^''^ 

For our problem, the model space consists of all of the parameters for the equation of state (EOS), /?,=i,...,a? , 
plus values for all of the masses of the neutron stars for which we have data, M;=i . at^^. From the pai^ameters. 
Pi, we can construct the EOS and solve the TOV equations to get a radius /?, for each of the neutron star 
masses M,. To be more concise in the following, we refer to our model M{p\,p2, ■■■, Pn,, Mi,M2, ..., M/v^^) 
as M(p\,„N,,,Mi,j^j^). Applied to this specific problem, eq. (|28] ) is 

P[Mipi...N,„M,,„M^)\D] = p[d\M(pi...n,,Mi,„n„)] 
X PlM{pL..N„^Mi,„N^)]nP[D\M]P[M]d''M\ (29) 

where N = Np + Nm is the dimensionality of our model space. The total number of EOS parameters is 
Np = S and the total number of neutron stars in our data set is Nm = 6. 

We construct our data D as a set of Nm probability distributions, Di{M, R) in the (M, R) plane, which 
are all normalized to unity, i.e. 

dM dRDi{M,R)=l'ii. (30) 

This normalization ensures that the data set for each neutron star is treated on an equal footing. We choose 
Mow = 0.8 Mq because current core-collapse supernova simulations fail to generate lower neutron star 
masses. The remaining limits, Mhigh = 2.5 Mq, /?iow = 5 km, and 7?high - 18 km are extreme enough to 
ensure that they have no impact on our final results. None of the probability distributions D inferred from 
the data described above has a significant probability for neutron st ars outside of these ra nges. Note also 



that the results for the neutron stai^s in Ml 3 and oj Cen are cutoff in I Webb & Barred (l2007|) for radii below 
8 km so our distributions also exhibit the same cutoff. 
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In order to apply equation ( |29l ) to our problem, we assume that P(D\JVl), the conditional probability of 
the data given the model, is proportional to the product over the probability distributions D, evaluated at the 
masses which are chosen in the model and evaluated at the radii which are determined from the model, i.e. 

p[d\M{pi...n,,Mi,„nJ]c^ Y\ Di{M,R)\M=M.,R=R(M.) (31) 

i=l...,NM 

This implicitly assumes that all of the data distributions D, are independent of each other and also indepen- 
dent of the model assumptions and prior distributions. Another required input for equation (1^91 ) is the prior 
distribution. We assume that the prior distribution is uniform in all of the np + hm model parameters, except 
for a few physical constraints on the parameter space described below. Taking a uniform distribution just 
means that the P{M) terms cancel from equation ( [291 ) and the integration limits become the corresponding 
prior parameter limits. 

In the maximum likelihood formalism, the function /'U)|At(pi...A?^^,Mi jv^^) is equivalent to the like- 
lihood function, and maximizing the likelihood function gives the best fit to the data. In the case that the 
pro bability distributions represen t Gaussian uncertainties, then the best fit is a equivalent to a least-squares 



fit (iBevington & Robinsonll2002h . In Bayesian inference, model parameters are determined using marginal 
estimation, where the posterior probability for a model parameter pj is given by 

P[Pjm{pj) = ^ fp[D\M] dpi dp2 ... dpj-i dpj^i . . . dpN^ dMi dM2 ... dMN^ (32) 

where V is the denominator in equation ( |29l ). without the model priors that determine the integration limits. 
The one-dimensional function P[pj\D]{pj) represents the probability that the j-th parameter will take a 
particular value given the observational data. Our problem thus boils down to computing integrals of the 
form in equation (l32l) . We use the Metropolis-Hastings algorithm to construct a Markov chain to simulate the 
distribution P u)|At(;?i..jVp»^i...WM) ■ For each point, we generate the EOS parameters p, and the neutron 
star masses M; from a uniform distribution within limits that are described below. The TOV equations are 
solved and this generates a R(M) curve and the radii for each neutron star, /?,. From these 6 masses and 
radii, the weight function P [D\M] is computed from eq. (|3T]) and the point is rejected or accepted using 
the Metropolis algorithm. In order to compute posterior probabilities P[pj\D]{pj), we construct several 
histograms of the integrand P[D|A1] from all of the points that were accepted. To construct the l-cr regions, 
we sort the histogram bins by decreasing probability, and select the first N bins which exhaust 68% of the 
total weight. A similar procedure is used for the l-cr regions. In order to constrain the full EOS of neutron 
star matter, we histogram the value of the pressure predicted by each EOS in the Markov chain on a fixed 
grid of energy density. To create the predicted curve, R{M), we construct a histogram of the predicted radii 
for each EOS in the Markov chain on a fixed mass grid. 

This analysis is easily extensible to a different number of EOS parameters or a different number of 
neutron star data sets. The only issue is that of computer time: the TOV equations must be solved for 
each point in the model space, and enough points must be selected to cover the model space fully. Another 
advantage is the explicit presence of the prior distributions, P{M). Although we have set these distributions 
to unity for this work, future work will utilize these terms to examine the impact of constraints on the 
equation of state from ten^estrial experiments. 
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4.2. Parameterization of the EOS 



We divide the EOS into four energy density regir nes. The region be low the transition energy density 
^trans * eo/2 is the crust, for which we use the EOS of iBaym et al.1 (119711) and iNegele & VautherinI (Il983l) . 
Here eo is the nuclear saturation energy density; it is convenient to remember that the nuclear saturation 
baryon density 0.16fm"^ corresponds to an energy density of k 160MeVfm"^^ and a mass density of ~ 
1.1 X 10^"* gcm~^. For etrans < e < £1, we use a schematic expression representing charge-neutral uniform 
baryonic matter in beta equilibrium that is compatible with laboratory data. Finally, two polytropic pressure- 
density relations are used in the regions e\ < e < 82 and e > £2- The densities ei and £2 are themselves 
parameters of the model. The schematic EOS for ejrans < e < ei is taken to be 

( K K' 3 1 

£^nBlmB + B+ —(u - if + (u - if + {I - Ixf [SkU^'"' + S pU^^ + -ficx{?>n^nbxfl^ \ (33) 

1 J- O J- \J^ t J 

where ns is the baryon number density, mg is the baryon mass, u = nslno, and x is the proton (electron) 
fraction. The saturation number density, uq, is fixed at 0.16 fm~^, the binding energy of saturated nuclear 
matter, B, is fixed at -16 MeV, and the kinetic part of the symmetry energy, Sk, is fixed at 17 MeV. The 
compressibility K, the skewness K' , the bulk symmetry energy parameter, S^ = S p + Sk (where Sp is the po- 
tential part of the symmetry energy), and the density dependence of the symmetry energy, y, are parameters, 
which we constrain to lie within the ranges specified in Table [3] These limits operate as constraints in our 
otherwise trivial prior distributions, P{M). To avoid bias in our results and to ensure that our model space is 
not over-constrained, we have intentionally made these ranges larger than normally expected from modern 
models of the EOS of uniform matter which ai^e fit to laboratory nuclei. While in principle the crust EOS 
for each set of EOS parameters could be different, as described in lSteinen (|2008r) . in practice the masses and 
radii are not strongly affected by changes in the crust at this level. The transition between the crust EOS and 
the low-density EOS is typically around half of the nuclear saturation density and is determined for each 
parameter set by ensuring that the energy is minimized as a function of the number density. We have opted, 
at this stage, not to include coiTclations between parameters that have been shown to exist from nucleai^ 
systematics or neutron matter calculations. For exampl e, the values of 5'v a nd y (or, equivalently, S^ and S s, 
the surface symmetry parameter) are highly correlated dSteiner et allllOOSh in liquid drop mass formula fits 
to nuclear masses. Such correlations will be considered in a future publication. 



The last term in equation (|33] ) is due to electrons. The proton fraction x is determined as a function of 
density by the condition of beta equilibrium 



de 



— = hcOn^nBxfl^ - 4 [Suu^^'^ + S pU^] {l-lx) = 0. 



dx 



(34) 



which has the solution 



where 



'=-A 



vl/3 



( Vt/ + 1 + 1) - ( VJ + 1 - 1) 



1/3 



d = 



2 f 



288 



he 



\Sku^i^ + Spuy\ 



(35) 



(36) 
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We also include muons in our equation of state, but they are not included in the above expressions for clarity. 
For lack of a clear theoretical understanding of the natu re of matte r abov e saturation densities, we 



parameterize the high-density EOS by a pair of polytropes. iRead et all (120091) have shown that such a 
parameterization can effectively model a wide variety of theoretical model predictions in this density range. 
Specifically, for energy densities above a parameterized value ei , the EOS is 

P = ^i£^+^^"' (37) 

where P is the pressure, n\ is the polytropic index, and ^i is a coefficient. Note that we are parameterizing 
the high-density EOS as polytropes in the total energy density e rather than the number density, and that 
this EOS is assumed to be in beta equiUbrium so it automatically includes leptonic contributions. Above a 
parameterized energy density, E2, we use a second poly trope 

P = K2£^^'^'"^ . (38) 

We choose the transition densities si and £2 and the polytropic indices ni and n2 as pai^ameters. The coeffi- 
cients Ki and K2 are determined by pressure and energy density continuity at the transition densities, with 
values hmited such that 1600 MeVfm"-' > £2 > ei > ISOMeVfm"^ Note that 1600MeVfm"^ is either 
larger than or nearly as large as the central energy density of most configurations. In addition, we limit 
£1 < 600 MeV fm~^ , so that the parameters of the schematic EOS maintain a close connection to their usual 
definitions in terms of properties near the saturation density. Finally, we limit the polytropic indices with 
0.2 < «i < 1.5 and 0.2 <n2 < 2.0; our results are insensitive to this choice of limits. 

While phase transitions are included in our models because they will appear- like successive polytropes 
with different indices, some parameter sets do not imply any phase transition, as they simply model equations 
of state which have effective polytrope indices which vary with density. Models with more than two strong 
phase transitions will not be well reproduced by our parameterization, but it is not clear that these models 
are particularly realistic. 

We thus have 8 total EOS parameters, K, K', Sy, y, ni, n2, £1, and £2, with limiting values summarized 
in Table [3l In addition to these parameter limits, some combination of parameters must be rejected because 
they ai"e unphysical. Unphysical combinations include those in which 

1. the maximum mass is smaller than 1.66 M©, which is 2-cr below the mass of PSR J 1903-1-0327, 1.74± 



0.04 Mo Champion et al.ll2008h : 



2. the EOS becomes acausal below the central density of the maximum mass star; 

3. the EOS is anywhere hydrodynamically unstable, i.e., has a pressure that decreases with increasing 
density; and 
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the maximum mass star has a maxim um stable rotation r ate less than 716 Hz, the spin frequency 
of the fastest known pulsar, Ter SAD (IHessels et alj|2006l) . T he spin frequency at which equatorial 
mass-shedding commences is given to within a few percent by (IHaensel et al.ll2009l) 



Ik - 1.08 



/ M 



1/2 



10 km 
R 



-ill 



kHz. 



(39) 



There has been claimed possible evidence for a higher spin frequency in XTE J1239-285 (|Kaaret et al. 



2007h . but this observation does not have strong statistical significance and has not been confirmed by 
subsequent observations. 



In addition to these criteria, during the Monte Carlo generation of neutron star mass sets for the Bayesian 
analysis, if one of the seven masses is larger than the maximum mass for the selected EOS, that realization 
is discarded and a new one is selected. 

To show that this model accurately represents significantly different EOSs, we fit it to the Skyrme 
model SLy4 which gives relatively small neutron star radii (of order 10 km for M = 1.4 Mq), and the 
field-theoretical model NL3 which gives rather larger radii, of order 15 km, for the same mass. These 
fits are shown in Fig. [5l illustrating that the associated EOSs are reproduced to within a few percent. Our 
parameterization includes many extreme models, including those like NL3 with rather large neutron star 
radii. We will find below that such models are ruled out, however, by the observational data. 



4.3. EOS Results from the Statistical Analysis 

We first analyze the EOS implied by observations for the case in which X-ray bursts are modeled as- 
suming Tph » R. The histograms for the EOS pai^ameters aie given in Figs. [6]-El The double-hatched 
(red) and single-hatched (green) regions outline the one- and two-sigma confidence regions which are sug- 
gested by the simulation. The one- and two-sigma parameter limits are summarized in Table |4] along with 
representative values obtained from terrestrial experiments. 

It is remai^kable that the inferred ranges for the schematic EOS pai^ameters K, K' , Sy, and y are consis- 
tent with the values derived from nuclear experiments. Especially significant is the inferred low value for y, 
a parameter which controls the pressure of the EOS in the range of 1- 3 eo and therefore the radii of neutron 
stars in the mass range 1 < M/Mq < 1.5 (ILattimer & PrakashI 12001 ). The value of y also controls the en- 



ergy and pressure of pure neutron matter. iHebeler & SchwenkI (|2009|) has shown that several recent studies 
indicate a convergence in predictions for neutron matter. For the saturation density used there, no = 0.16 
fm~^, the mean values of the neutron matter energy and pressure are 16.3 MeV and 2.5 MeV/fm^^, which 
imply 5,, = 32 MeV and y = 0.28 with the parametrization of Eq. (1331 ). well within our predicted ranges. 



One must take care in comparing experimental constraints directly to our results, however. The param- 
eters determined experimentally are often "local" quantitie s, in the sense that t hey are properties of the EOS 
only at densities close to saturation. Also, the results from lTsang et al.l (l2009|) mostly constrain the EOS in 
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Table 3. Prior limits for the EOS parameters 



Quantity 


Lower limit 


Upper limit 


K (MeV) 


180 


280 


K' (MeV) 


-1000 


-200 


Sv (MeV) 


28 


38 


r 


0.2 


1.2 


«i (fm-3) 


0.2 


1.5 


«2 (fm-3) 


0.2 


2.0 


ei (MeV fm"^) 


150 


600 


£2 (MeV fm^-^) 


£i 


1600 



10^ 


^__.^;;::::::=;;s=ss= 




E NLs,,..--^;^^::!;^^^^^ 




- /^ XSLy4 


CLh 


- // 


10 


i 


1 






T 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
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8 (fm"^) 



Fig. 5. — Fits of our parametrized EOS (dotted lines) to the Skyrme EOS SLy4 and the field-theoretical EOS 
NL3{solid lines). 
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region near half the saturation density. We have used the schematic EOS parameters over a somewhat larger 
range of densities, from 1/2 to, typically, up to 2 or 3 times the saturation density. 

The predicted skewness for the schematic equation of state has a relatively small magnitude, centered 
at K' = -280 MeV. The most probable value for the transition between the schematic EOS and the first 
polytrope is around twice the saturation density, sq. The index of the first polytrope is sharply peaked 
around 0.5, which corresponds to a polytropic exponent yi = I + l/«i ^ 3 which is rather large. Coupled 
with the small magnitude of the skewness, this implies a quite stiff EOS at supernuclear densities. Finally, 
note that, typically, nz > n\, indicating that the EOS softens at high densities. A summary of the results of 
the EOS parameters for both the schematic EOS and the poly tropes are given in Table ID 

We next consider the results of parameter fitting when the X-ray burst data is modeled assuming 
Tph - R. Results are summarized and compared to the previous case in Table |4] Although we don't show 
histograms for the parameters in this case, their behavior is similar to those of Fig. [6] and |7] modulo the 
different means and variations of the two cases. It is significant that the small value for y found for the case 
Tph » /? is duplicated for rph = R, supporting a weak density dependence for the sym metry energy and a 



consequent small estimate for neutron star radii in the mass range 1-1.5 Mq. We note that lOzel et al.l (120 lOT) 



also conclud ed, on the basis o f observations of X-ray bursts, that the implied neutron star radii were rela- 



tively small. lOzel et al.l (|2010l) obtained radii approximately 1 km smaller than we do, however, for rph = R 
(which are already 0.5-1 km smaller than those we obtain for rph » R) because their analysis favors the 
high-redshift solution, and because they do not impose the same causality constraint on their Monte Carlo 
sampling. 

Although the low-density EOS is not strongly affected by the assumption that rph - R, at higher 
densities this choice leads to a softer EOS: the magnitude of the skewness K' and first polytropic index n\ 
are both larger for the case rph - R (Table |4l). The differences between these predicted EOS can be easily 
seen by referring to Fig. [8] Each panel displays an ensemble of one-dimensional histograms over a fixed 
grid in one of the axes (note that this is not quite the same as a two-dimensional histogram). The bottom 
panels present the ensemble of histograms of the pressure for each energy density. This was computed in 
the following way: for each energy density, we determined the histogram bins of pressure which enclose 
68% and 95% of the total MC weight. The locations of those regions for each 1-D histogram are outlined 
by dotted and solid lines, respectively, and these form the contour lines. These 1- and 2-cr contour lines give 
constraints on the pressure of the EOS as a function of the energy density as implied by the 6 neutron star 
data sets and are presented in Tables |5] and [6] along with the most probable EOS for the cases rph » R and 
rph - R, respectively. The softer nature of the EOS in the rph - R case is most apparent at the highest energy 
density, for which the pressure is 10% less than in the rph » R case. 

The upper panels of Fig. [8]display the ensemble of histograms of the ratio of the pressure of the EOS to 
the pressure of SLy4 over a grid of energy density. The choice of the SLy4 EOS here is essentially arbitrary, 
and just assists in plotting the results since the pressure varies over a couple orders of a magnitude. The 
infeiTcd pressure ratio from several Skyrme models is also shown. These Skyrme models were selected in 



Steiner & WattsI (120090 because they are a representative set of the recommended models from IStone et al. 
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Fig. 6. — Histograms for the compressibility K, skewness K', symmetry energy Sy,, and density dependence 
of the symmetry energy y. The l-cr (double-hatched) and l-cr (single-hatched) confidence regions are also 
indicated. These results assume rph » R. 
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Fig. 8. — Probability distributions from the Bayesian analysis of the the parametrized EOS. Upper panels 
show the ratio of the most probable EOS to the pressure of SLy4. The lower panels show the pressure of the 
most probable EOS. The left panels show results under the assumption rph = R, and the right panes show 
results assuming rph » R. In all panels, the solid (dotted) contour lines in each panel show the 2-cr (1-cr) 
contours implied by the data. 



■28- 



(120031) which have symmetry energies which do not become negative at high density. In both the rph = R 
and Tph » R cases, there appears to be a softening of the EOS at low densities which is incompatible with 
some of the Skyrme models. Some models are apparently ruled out independ ently of assumptions a bout the 



radius of the photosphere in X-ray bursts: field-theoretical models hke NL 3 (ILalazissis et al.ll 19971) . which 



have stiff" symmetry energies; FSUGold (ITodd-Rutel & Piekarewica 120051) . which has a softer symmetry 



energy but becomes too soft at high density; and APR (jAkmal et al.lll998h . which becomes too stiff at high 
densities. The neai^ly vertical nature of APR for energy densities just over 200 MeV fm~^ is due to the phase 
transition in APR to matter which includes a n^ condensate. Even after applying a Gibbs phase construction, 
the pressure increases only very slowly with density in the small mixed phase region. 



4.4. Mass and Radius Results from the Statistical Analysis 

The upper panels in Figure |9]present our results for the predicted mass-radius relation according to our 
two assumptions regarding the photospheric radius for X-ray bursts. They give the ensemble of histograms 
of the radius over a fixed grid in neutron star mass. The width of the contours at masses below 1 M© tends to 
be large because the available neutron star mass and radius data generally implies larger masses. In general, 
the implied M-R curve suggests relatively small radii, which is consistent with our conclusions regarding 
the softness of the nuclear symmetry energy in the vicinity of the nucleai^ saturation density. Tables |7] and [8] 
summarize the 1- and l-cr contour lines from these panels, and give as well the most probable M-R curve. 
The assumption that r^h - R implies smaller radii for neutron star masses greater than about 1.3 M© and 
perhaps very small radii for masses in excess of 1.5 Mq. This is suggestive of the onset of a phase transition 
above the nuclear saturation density or perhaps simply the approach to the neutron star maximum mass limit 
in this case. The mass and radius contour lines, determined in the same way as for Figure [H are also given 
in the figure, and summarized in Tables |7] and [H The radii in the rph = R case average about 1 km less than 
in the rph » R case, except around 1.4 Mq, where the difference is about 0.4 km. 

The lower panels in Figure |9] summarize the output probability distributions for M and R for the 6 
neutron stars which were used in this analysis. Note that the scales have been modified, and these output 
probability distributions are much smaller than the input distributions given in Figures [H [2l and HI The 
combination of several neutron star mass and radius measurements with the assumption that all neutron 
stars must lie on the same mass-radius curve puts a significant constraint on the mass and radius of each 
object. Note that these output probability distributions closely match the implied M vs. R curves presented 
in the upper panels of Figure|9l The tendency for smaller radii when M > 1. 3-1.4 M© is apparent. 

The M and R constraints for each object are given in Table|9l with their corresponding l-cr uncertainties. 
The three PRE burst sources suggest masses near 1.5 Mq and the radii near 1 1 km in the case where rph = R. 
We have already observed that this agreement is due to the extreme restrictions placed on the acceptability 
of points during the Monte Carlo sampling. In the case rph » R, the PRE burst sources are predicted to have 
a wider range of masses, from 1.3 to 1.6 solar masses. The quiescent LMXB masses tend to be smaller, 
but are strongly dependent on assumptions about the radius of the photosphere in the PRE burst sources. 
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Fig. 9. — The upper panels give the probabihty distributions for the mass versus radius curves implied by 
the data, and the solid (dotted) contour Unes show the 2-o" (1-cr) contours implied by the data. The lower 
panes summarize the 2-cr probability distributions for the 6 objects considered in the analysis. The left 
panels show results under the assumption rph = R, and the right panes show results assuming rph » R. The 
dashed line in the upper left is the limit from causality. The dotted curve in the lower right of each panel 
represents the mass-shedding limit for neutron stars rotating at 716 Hz. 



30 



Particularly uncertain is the mass for X7. The observations imply a rather large value of Rco, which is 
compatible with a small radius if the mass is large (rph » R), or a large radius if the mass is small (rph = R). 
The stars in Ml 3 and co Cen show the opposite trend. They have small values of Rco which is compatible 
with small radii in the rph » R case if the mass is relatively small. In general, the predicted radii of all stars 
range from about 10 km to 12.5 km, with the exception of X7 which may have a large radius. Note that a 13 
km radius for an 0.8 M© star is beyond the limit implied by rotation at 716 Hz and thus if the neutron star in 
X7 is observed to spin rapidly enough, a much larger mass and smaller radius would be implied instead. 

The largest difference between the predicted equations of state between the rph = R and rph » R 
cases is the high-density behavior. This leads to large differences in the predicted maximum masses. The 
probability distributions for the maximum neutron star mass are given in Figure[TOl along with the associated 
1-cr confidence regions. The two probability distributions are arbitrarily normalized so that their peak is 
unity. These results are strongly dependent on assumptions of the photospheric radius at touchdown. The 
two-peaked behavior in the case rph = R suggests a possible phase transition could match the data and 
im plies a max i mum mass very close to the observed limit of 1.66 M©. This result is similar to that claimed 
by lOzel et al.l (|2010l) . but there it is stated that the results are incompatible with a nucleonic equation of 
state. Our results do not support this extreme interpretation. Although the neutron star radii implied by this 
analysis ai^e small, they are not small enough to require a phase transition; for neu tron stars of mass 1 .4Mo, 
radii smaller than 10 km can be generated by purely nucleonic equations of state dSteiner et al.ll2005r) . 



5. Discussion 



In this paper, we have determined an empirical dense matter equation of state from a heteroge neous 
datas e t containing PRE bursts and quiesce n t thermal ernission f rom X-ray transients. Previous works (|Qzel 
20061 : bzeletal-lboogl : lOuver et allboiol : loiiver et allboiObh have demonstrated the potential utiUty of 
PRE bursts for determining neutron star masses and radii. Their analysis assumes that the photosphere 
at touchdown has fully retreated to the quiescent radius, and that the high-redshift solution is favored. We 
argue that their model requires re-examination. First, we find no grounds on which the high-redshift solution 
in the case rph = R could be favored. Second, internal consistency (i.e., the obUgatory rejection of an 
overwhelming number of Monte Carlo realizations) implies the assumption that rph = R is suspect and 
should be generalized. We then explored an alternative, that the radius of the photosphere is extended and 
does not recede until later in the observed burst. A larger photosphere indeed provides internal consistency 
without requiring strong cuts on the observed values for flux, normalization, distance, and composition. 



There are other sources of systematic enors in the PRE burst model. iBoutloukos et al.l (120101) found 
that RXTE/PCA spectra of bursts from 4U 1 820-30 were better fit by Planck or Bose-Einstein spectra rather 
than Comptonized spectra with large color correction factors. If the color correction factor is indeed of order 
unity, then the high color tem peratures would indicate a locally super-Eddington flux, even in the tail of the 
burst (|Ebisuzaki et al.lll984) . In addition, at near-Eddington fluxes there are other possible complications. 
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Table 4. Most probable values of the EOS parameters and their associated l-cr uncertainties. 



Quantity rph » 7? rph = /? Experiment 



K (MeV) 
-K' (MeV) 
S (MeV) 

r 



Schematic EOS parameters 



216:^3 

280!4io 

90+5.4 

26+°i5 



190+^2 
500+290 

-^" -170 

26+"-^'' 



230-250 (Colo et al. 2004) 

28-34 (Klupfel et al. 2009) 
0.35-1.0 (Tsang et al. 2009) 



High-density EOS parameters 



ni 

ei (MeV/fm^) 
£2 (MeV/fm3) 



Qcj+0.12 

1 90+0.49 

^•^-'-0.16 

290+80 

"^ -170 



67+°-20 
"■" -0.20 

1 90+0.59 

-^-0.22 
320+^20 

880!«o 






U 0.5- 




Fig. 10. — The probability distributions for the maximum mass. Because of the observation of a neutron 
star of at least 1.66 M©, we reject all curves for which M^ax < 1.66 M© and thus the probability is cutoff 
below this value. The shaded regions indicate the l-cr confidence regions. 
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Table 5. Most probable values and 68 and 95% confidence limits for the pressure as a function of energy 

density, for rph = R. 



Energy density 2-o" lower limit 1-cr lower limit Most probable 
MeV/fm^ MeV/fm^ 



1-cr upper limit 2-cr upper limit 



150 


1.97 


2.02 


2.16 


3.77 


4.11 


200 


4.57 


5.01 


5.25 


7.37 


8.26 


250 


8.03 


9.24 


11.99 


12.63 


13.89 


300 


11.68 


14.23 


17.48 


18.88 


20.76 


350 


16.23 


20.00 


22.47 


25.95 


28.72 


400 


22.73 


26.72 


31.96 


33.82 


37.29 


450 


30.51 


34.73 


40.75 


43.69 


48.43 


500 


38.28 


43.97 


49.72 


56.10 


63.79 


550 


47.44 


53.93 


61.06 


71.02 


84.56 


600 


58.92 


65.10 


72.64 


87.99 


110.8 


650 


70.68 


77.15 


82.31 


107.2 


139.4 


700 


84.24 


90.40 


97.50 


128.9 


172.2 


750 


99.90 


105.2 


117.6 


153.4 


207.4 


800 


115.6 


121.0 


132.3 


179.6 


239.2 


850 


132.5 


136.5 


155.8 


206.1 


270.3 


900 


151.0 


155.8 


173.8 


263.0 


304.9 


950 


169.7 


177.5 


191.6 


302.2 


338.5 


1000 


190.8 


197.9 


207.7 


337.3 


371.1 


1050 


213.2 


220.1 


236.5 


364.2 


406.2 


1100 


234.8 


243.0 


253.5 


389.8 


439.3 


1140 


257.2 


268.1 


277.5 


423.2 


476.7 


1200 


280.2 


294.5 


308.3 


459.7 


516.3 


1250 


300.5 


321.1 


343.4 


462.8 


551.2 


1300 


323.4 


346.2 


379.2 


504.2 


585.8 


1350 


352.0 


369.9 


411.0 


528.2 


622.8 


1400 


381.2 


398.4 


436.4 


561.5 


666.9 


1450 


402.4 


432.6 


453.1 


612.5 


712.3 


1500 


427.1 


460.4 


502.1 


643.9 


749.2 


1550 


459.3 


484.6 


518.8 


676.0 


799.6 


1600 


480.0 


520.7 


550.8 


718.2 


848.6 
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Table 6. Most probable values and 68 and 95% confidence limits for the pressure as a function of energy 

density, for rph » R. 



Energy density 2-o" lower limit 1-cr lower limit Most probable 
MeV/fm^ MeV/fm^ 



1-cr upper limit 2-cr upper limit 



150 


1.84 


2.00 


2.15 


2.52 


3.25 


200 


4.44 


4.97 


5.76 


6.06 


6.88 


250 


8.10 


9.14 


10.67 


11.50 


12.71 


300 


12.34 


14.61 


15.87 


19.08 


21.43 


350 


18.19 


21.52 


24.43 


29.16 


34.17 


400 


25.32 


30.46 


35.80 


42.86 


52.07 


450 


34.56 


41.54 


49.11 


60.99 


77.97 


500 


45.32 


55.19 


69.75 


84.54 


107.7 


550 


57.47 


73.14 


89.83 


114.1 


138.1 


600 


70.62 


94.02 


115.5 


145.0 


168.1 


650 


85.07 


116.8 


149.7 


175.8 


199.8 


700 


100.8 


141.0 


173.2 


206.8 


230.6 


750 


115.7 


165.8 


203.3 


237.8 


263.8 


800 


131.5 


190.6 


232.7 


268.8 


297.7 


850 


154.0 


214.9 


254.1 


299.9 


330.9 


900 


175.0 


239.7 


277.4 


331.6 


364.7 


950 


196.7 


264.4 


310.0 


364.2 


399.7 


1000 


217.0 


289.9 


336.2 


397.8 


437.6 


1050 


236.5 


316.5 


371.9 


432.7 


474.1 


1100 


260.3 


343.5 


414.4 


468.9 


509.8 


1140 


282.1 


369.6 


458.2 


504.7 


547.6 


1200 


307.3 


395.4 


457.7 


541.2 


592.7 


1250 


333.2 


424.7 


494.3 


581.2 


639.6 


1300 


356.8 


454.5 


550.4 


622.4 


680.9 


1350 


381.4 


481.1 


603.9 


660.6 


726.2 


1400 


407.7 


512.5 


608.5 


704.7 


777.2 


1450 


434.2 


542.8 


672.0 


749.1 


824.9 


1500 


460.5 


572.1 


661.5 


791.8 


877.5 


1550 


487.5 


606.4 


732.5 


840.2 


932.8 


1600 


514.9 


635.2 


734.9 


884.5 


979.8 
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Table 7. Most probable values and 68 and 95% confidence limits for neutron star radii of fixed mass, for 

Mass 2-cr lower limit 1-cr lower limit Most probable 1 -cr upper limit 2-cr upper limit 
Mq km 



1.0 


11.06 


11.35 


12.09 


12.62 


13.31 


1.1 


10.97 


11.28 


11.87 


12.35 


13.00 


1.2 


10.90 


11.23 


11.69 


12.11 


12.65 


1.3 


10.81 


11.13 


11.42 


11.85 


12.28 


1.4 


10.68 


10.98 


11.32 


11.60 


11.93 


1.5 


10.42 


10.75 


11.04 


11.37 


11.65 


1.6 


9.93 


10.37 


10.87 


11.17 


11.47 


1.7 


9.44 


10.05 


10.65 


11.12 


11.42 


1.8 


9.60 


10.14 


10.65 


11.06 


11.38 



Table 8. Most probable values and 68 and 95% confidence limits for neutron star radii of fixed mass, for 

rph » R. 

Mass 2-cr lower limit 1-cr lower limit Most probable l-cr upper limit 2-cr upper limit 
Mq km 



1.0 


11.29 


11.75 


12.09 


12.30 


12.57 


1.1 


11.31 


11.71 


11.94 


12.28 


12.51 


1.2 


11.30 


11.65 


11.96 


12.24 


12.47 


1.3 


11.24 


11.58 


11.81 


12.21 


12.45 


1.4 


11.13 


11.49 


11.83 


12.18 


12.45 


1.5 


10.94 


11.39 


11.83 


12.17 


12.45 


1.6 


10.63 


11.30 


11.70 


12.17 


12.49 


1.7 


10.42 


11.21 


11.70 


12.13 


12.54 


1.8 


10.43 


11.10 


11.57 


12.05 


12.47 
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such as a photon bubble instabiUty (IHsu et all 19971) that may impact the observed properties of a PRE X-ray 
burst. 



A recent analysis of PRE bursts in 4U 1724-307 by ISuleimanov et al.1 (120101) also found that the pho- 
tosphere is somewhat extended at touchdown, rph w 27?. Using a new model of the atmosphere (the details 
of which are not yet published), they obtain masses and radii that are both larger than we would have found. 
This result is principally due to a larger value of the color correction factor, f c, predicted by thei r model 
atmosphere. It is 10-15% larger than what one would obtain from the model of iMadej et al.l (120041) for the 



same values of f /f Edd> surface gravity g, and composition. In thi s case, we find values of a oz f^^ which 



are about 25% larger and /?oo «: f^ that are about 25% smaller than ISuleimanov et all (120101) did. Assuming 
Tph » R, the predicted mass scales with f^. When we analyse 4U 1724-307 with the lower value oi fc, we 
get masses and radii that are consistent with those from the other three PRE burst sources. This highlights an 
important avenue for future work, namely that a clearer understanding of the atmosphere and, in particular, 
fc are essential. 

Further progress in using PRE bursts to constrain neutron star masses and radii clearly requires better 
models of the spectral evolution during X-ray bursts. It is important to note that increases in the precision of 
mass and radius estimates are not absolutely necessary, as we have shown that existing errors, large as they 
ai^e, do not inhibit placing interesting limits to the equation of state when combined with other observations. 
It is very important, however, to resolve the uncertainty regarding the location of the photospheric radius 
at "touchdown" and to characterize systematic uncertainties, including potential correlations between Fjd. 
A, fc, and X in the spectral models. Also, should a larger range in fc be required, then the uncertainties on 
masses and radii will increase accordingly. 

Our results imply that the EOS in the vicinity of the nuclear saturation density is relatively soft. As 
shown by lLattimer & PrakashI (120011) . this is primarily due to the weak dependence of the nuclear symmetry 
energy with density. This conclusion is robust with respect to variations in how PRE burst sources are 
modeled, and is strengthened by the small values of /?oo deduced for the globular- cluster LMXB sources in 
M13 and cj Cen. This conclusion has immediate and important ramifications for laboratory nucleai^ physics, 
in particu lar for the scheduled parity-violating electro n scattering measurement of the neutron skin thickness 
of 208pb ^Michaels et allbood : Horowitz et allboOlb . In the contex t of the liquid dropl et model, the ratio of 
the surface and volume symmetry coefficients can be expressed as dSteiner et al.ll2005h 






9£.so h 



(40) 



where £"^0 - 19 MeV is the surface energy coefficient for symmetric nuclei. The integrals I\ and I2 for the 
schematic EOS described by equation (l33l) ai^e given by 



Jo 

I2 = I ^/u 

Jo 



(1 - u) du. 



Sku^l^ + {S,- SkW 



1 



1 



1 - u VI - a(l - u) 



du. 



(41) 
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where a - K' /(9K). The predicted neutron excess in the center of a nucleus {N,Z) is then 

-1 



N-Z 
N + Z 



1 + 



Ss 



5,,Ai/3 



and the neutron skin thickness is 



3 2roS, 
SR= J-- 



(42) 



(43) 



5 3 5vl-52' 

where r^ - (47rno/3)~'''^. With a value of y = 0.26 and a - -0.141 and errors as established in Table 4, we 
deduce S JSy ^ 1 .5 ± 0.3 and 6R(^^^Pb) - 0. 15 ± 0.02 fm, a value at the lower end of expectations. 

We also conclude that in our preferred model, in which the photospheric radius is not restricted to be the 
stellar radius, the EOS at high densities is stiff enough to support a maximum mass of order 2 Mq or greater. 
This result is supported by the simultaneous observations of LMXB sources with both small and large values 
of Rco, if one rejects the possibility of the existence of neutron stars with masses less than about 1 Mq. The 
latter condition seems to be borne out by models of massive stai^ supernovae, which strongly suggest that 
proto-neutron stars are born with high trapped lepton fractions and moderate entropies. Hydros tatic stability 
requi res, in this case, that protoneutron stars are bound only if their mass exceeds about 1 M© (IStrobel et al. 
19991) . ruUng out the existence of cold, catalyzed neutron stars of smaller masses. Had this assumption been 
used as a prior condition in our analysis, the LMXB sources with small values of Ra, would further support 
solutions with EOS parameters favoring large maximum masses. 

Obviously, more neutron star observations with mass and radius constraints would enable one to im- 
prove our constraints. It is a particular advantage of our methodology that new observations can be integrated 
into the formalism easily, i f estimates of a two- dimensional probability distribution of mass and radius can 
be made. As pointed out by lHeinke et al.l (120061 ). it is particularly important that atmosphere models treat the 
surface gravity self-consist ently. In particular, spectral features would be very constraining, and although 
the rn ost recent result from ICottam et al.l (l2002h has not been verified in longer obser vations dCottani et al. 
2008h . future determinations of the surface gravity may provide strong constraints. IHo & Heinkd (l2009h 
also demonstrated that constraints on the mass and radius of the Cas A supernova remnant may be obtained. 
These constraints are roughly consistent with our rph » R results described above. 

Although we did not include the estimate of the mass and radi us (M = 1.7 ± 0.3 M© and R = 
11.5 ± 1.2 km) of RX Jl 856-3754 (jPons et al.ll2002l : Iwalter et allboioh in our baseline fit, we found them 
to be remarkably consistent with those of the other 6 neutron stars. All of the results we obtained for the 
parameters of the EOS, the pressure versus energy density curve, the mass versus radius curve, our estimate 
of the maximum mass, and the predicted individual masses and radii for all 6 other sources were unchanged 
to within l-cr. For example, the most probable radius of a 1.4 M© star would change by less than 0.01 km 
and the most probable pressure at an energy density of lOOOMeVfm"^ would increase by 0.03%. 

In addition to X-ray bursts, quiescent low-mass X-ray binaries, and thermally-emitting isolated neu- 
tron stars, other potentia l methods for determining neut r on stai' masses and ra dii exist. These include 
neutron star seismology (ISamuelsson & AnderssonI l2007l : ISteiner & WattsI 120091) . pulse profiles in X-ray 
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pulsars (ILeahv et al. 



2009; 



Morsink & Leahyll2009l). and moment of inertia measurements in relativistic 



binaries (ILattimer & Schuta 120051) . Grav itational wave signals of neutron star mergers may also provide 



significant constraints (IFerrari et al.ll2010l) . Measurements of the thickness of the neutron star crust, which 
controls the crust cooling of transiently accreting LMXBs dShtemin et al.ll2007uBrown & Cumniindl2009l) 



as we ll as the distribution of observed cooling curves compared to a minimal cooling model (e.g. JPage et al 
2009h . could also constrain neutron star masses and radii. 



The Bayesian analysis in Section |4] is a novel procedure for combining mass and radius constraints 
from disparate objects to form new constraints on the mass versus radius curve and the EOS of d ense matter. 



The construction of the EOS f rom astrophysical obse r vations has also be en recently addressed in Read et al. 



(l2009h . lOzel & PsaltisI (120091) . and lOzelet all (120101) . Read et all (120091) showed that piecewise polytropes 
with relatively few parameters could effectively describe more sophisticated models for a high-density EOS, 
but confined attention to constraints stemming from observations limiting the neutron star maximum mass 
and maximum spin rate. They did not attempt reconstruction of the EOS from simultaneous mass and 
radius measurements. We therefore obtain stronger cons traints on the EOS, bu t have utilized observational 
data which contains significant systematic uncertainties. lOzel & PsaltisI (|2009h examined the constraints on 
the EOS obtained from a synthetic data set obtained with simultaneous mass and radius measurements of 
either 5% or 10% accuracy for thi^ee separate objects; they found that although individual EOS parameters 
were difficult to estima te precisely, significan t congelations between them could be obtained. The Jacobian 
technique employed in lOzel & PsaltisI (120091) is a simple approximation of our full Markov Chain Monte 
Carlo method. That Jacobian technique requires that the number of EOS parameters and the number of 
neutron stars are equal, and thus provides only an incomplete marginalization over the EOS parameters. 

One could imagine various alternatives to the statistical procedure in Section |4| One possibility is 
the use of prior distributions to represent either astrophysical or nuclear input. Neutron star masses in 
some double-pulsar systems are well-measured, and could provide a prior mass distribution. It is not clear, 
however, that either low-mass X-ray binaries or isolated neutron stars follow this same initial mass function. 
As mentioned above, there are strong theoretical reasons to suspect that neutron stars cannot be formed with 
less than about 1 .0 Mq. In the analysis above, we chose 0.8 M© to ensure that the boundary does not interfere 
with results near the proto-neutron star minimum mass. There are several nuclear physics observables 
which can be used to constrain the EOS through measurements of properties like the compressibility and 
the symmetry energy. Several authors have also computed the EOS of neutron matter directly from nucleon- 
nucleon interactions (|To1os et al.ll2008l: iHebeler & SchwenkI 120091) . and we have found that our predicted 
EOS is consistent with these studies. In the context of this work, information from calculations of the neutron 
matter EOS, neutron ski n thicknesses, the s urface symmetry- volume symmetry energy correlation observed 
from mass formula fits dSteiner et al.ll2005h . information from giant dipole resonances, and so forth, could 
also provide constraints for the schematic EOS parameters described above. We have chosen not to use this 
information in this work, in part to ensure that our results are free from extra model dependencies. Future 
work on implementing more nuclear physics input into the analysis of the observational data is certainly 
warranted. 
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One might also question, for the high-density EOS, whether or not our prior distributions for the poly- 
tropic indicies ought to be treated as uniform, as a uniform distribution in the polytropic index is different 
than, for example, a uniform distribution in the pressure at any fixed energy density, or in the polytropic 
exponent. Finally, one could consider reformulating the model directly in terms of the observables like flux 
and distance instead of the "two-step" procedure we have used above which uses a MC simulation to gen- 
erate masses and radii from the observables for each neutron star; and then afterwards performs a Bayesian 
analysis from the output of these initial MC results. This alternative would strongly disfavor rph = R, since 
so many MC realizations must be rejected in this picture. 

It is important to note that the constraints we have obtained are a guide, but will likely require revision 
in the future. Because we only used 6 mass and radius measurements a single additional measurement could 
have a significant impact on our results. Alternatively, if one of the 6 input probability distributions used 
above changes significantly, because, for example, of a new understanding of systematic uncertainties, our 
final constraints on the EOS would change accordingly. We have already noted that there is significant 
tension on our results from assumptions about the photospheric radius of X-ray bursts, and also from the 
fact that the neutron star in Ml 3 has a small value of /?« while X7 in Terzan 5 has a much larger value of 

It is possible that there exist extreme models which live in very small regions of parameter space and are 
not fully sampled in this work. One example of such an EOS would be those which exhibit the "twinning" 
phenomenon, i.e. the p resence of a turning point in the mass vs. radius curve which admits a new family of 



compact neutron stars (|Glendenning & Kettnenl2000l : ISchaffner-Bielich et al.ll2002h Such solutions should 



be addressed in future work. 



Ozel et al.l (120101) claim that constraints from the PRE burst sources imply that equations of state which 
contain only nucleonic degrees of freedom are inconsistent with data. We disagree with their conclusion. 
Although we find the EOS must be quite soft at moderate densities, we do find many nucleonic models 
which are consistent with the data. More importantly, however, we find that the conclusion of extreme 
softening evaporates if we use a slightly different model for the PRE burst sources that accounts for sys- 
tematic uncertainties. The inclusion of mass and radius data from other neutron star sources supports our 
interpretation of the PRE burst sources. While our results do not rule out a phase transition at supernuclear 
densities, extreme softening of the EOS is not compatible with observations. Rather, the implication is that 
the maximum mass is likely large, greater than 1.8 M©. 
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A. Observations of Type I X-Ray Bursts 

In this section, we describe in detail our probability distributions for the angular emitting area, or 
normalization, the distance, the touchdown flux, and the photospheric hydrogen mass fraction for the three 
sources with PRE bursts used in this analysis. We also note, where appropriate, how the distances compare 
to the value obtained by assuming the maximum flux is less than the Eddington value. 



z,<,^,.n.ikpc 

t f^^ max J 



M \/10 ^ ergs cm ^s '^"' 



1.4 Mp 
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A.l. EXO 1745-248 



Nonnalization: The normalization factors obtained from the two PRE bursts in lOzel et all (l2009h are 
A/(l km^ kpc~^) - 1.04±0.01 and 1.30±0.01. These two measurements differ significantly, and this means 
that either the color correction factors for the two bursts differ by roughly 6% o r the geometry of the two 



bursts is different. It is not clear how to resolve this conflict. lOzel et al.l (120091) use a boxcar distribution 
with A = 1.16km^kpc"^ and AA = 0.13 km^kpc"^, which represents the choice with the minimum possible 
uncertainty. We instead choose a Gaussian centered at A = l.lVkm^kpc"^ with cta = O.lBkm^kpc"^, 
which ensures that the two observations are included to within 1 cr. This choice is consistent with the 
objects below, which also have normalization factors which are simulated with Gaussian distributions, but 
with notably smaller uncertainties. 



Distance: EXO 1745-248 is located in the globular cluster Terzan 5. lOrtolani et all (120071) ana- 
lyzed the distance to Terza n 5 using both the NICMOS instrumental magnitudes and the calibrations from 
Stephens et all (120001) and ICohn et all (l2002h . The combined distance estimation is D = 5.9 ± 0.9 kpc, 
where the uncertainty has been obtained from the standard deviation of the three different distance measure- 
ments, suggesting a Gaussian with the same standard deviation. On the basis that the distance measurement 
from the NICMOS i nstrum ental magnitudes was independent of photometric calibrations and thus more 



accurate, lOzel et all (|2009r) used the NICMOS distance of D = 6.3 kpc with AD ^ 0.32 kp c. This choice 



assigns, however, a zero probability to a distance of 5.9 kpc, the central value suggested in lOrtolani et al. 



(l2007h . Until this distance measurement is more clearly determined, we choose a Gaussian distribution with 
D = 6.3 kpc and ctd = 0.6 kpc. 



Touchdown flux: We employ the result from lOzel et al.l (|2009h . which is a Gaussian distribution with 



-2 ,-2 



^TD,oo = 6.25 ergs cm s and o-p = 0.2 ergs cm s 



Hydrogen mass fraction: iGalloway et all (l2008h noted that EXO 1745-248 has exhibited long Type 
I bursts with estimated values of jdt ^persistent/ jdt f burst * GM/R/Enac ~ 20-46. These long bursts 
do not always show a strong thermal evolution (Galloway, private communication); but if they are indeed 
thermonuclear in origin, then the low ratio of persistent to burst fluence indicates H-rich fuel and larger 
values of X. In contrast to the determination of X from these long bursts, this object has also been identified 
as an ultracompact binary in iHeinke et al.l (|2003h . through a phenomenological assessment of its spectral 
properties, suggesting that X is small, < X < 0.1. The radius expansion bursts have short durations, which 
would be consistent for ignition of a pure He layer (i.e., the hydrogen is consumed by the stable hot CNO 
cycle). We retain the full range of hydrogen composition, < X < 0.7 . 



A.2. 4U 1608-522 



Normalization: The normalization factors f or the 4 PRE burst s are A /(I km^ kpc ^) = 3. 267 ± 0.047 



3.302 ± 0.049, 3.258 ± 0.054, and 3.170 ± 0.047 (iGuver et al.ll2010t) . We use the average from iGuver et all 



torn . 3.246 ± 0.024. 
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Distance: iGiiver et all (120101) give a Gaussian distribution with D = 5.8 kpc and ctd = 2.0 kpc with a 
cutoff below 3.9 kpc. We also employ this result. We note that if the flux is indeed less than the Eddington 
value, then D < 4.36kpc(M/1.4Mo)'^^ for the central value of the touchdown flux. 

Touchdown flux: Two of the four PRE bursts gave a value for t he touchdown flux, f td.oo = (15.58 + 
0.82) ergs cm"^ s^^ and (15. 14±1.05)ergscm~^ s"^ (jouver et alboioh . The value (15.41±0.65) ergs cm"^ s"^ 
was obtained from the fit in 



GuveretaLlteOld) . 



Hydrogen mass fracti on: The bursts in 4U 1608-522 suggest an accretion rate in of 3%-5% MEdd, 
which suggests H ignition (|Galloway et al.ll2008h : the brighter bursts from this system are of short duration, 
however, so it is likely that much of the hydrogen is consumed via the HCNO cycle prior to He ignition. As 
with EXO 1745-248, we use the fuU range < X < 0.7. 



A.3. 4U 1820-30 

Normalization: The three PRE bursts for which a n ormalization was ob tained give A/(l km^ kpc~^) = 
0.8886 ± 0.0373, 9668 ± 0.0339 and 0.9040 ± 0.0200 (JGuver et allboiObh . We use the value quoted in 
buyer et al.l(l2010bh . 0.9198 ± 0.0186. 



Distance: 4U 1820-30 is in the globular cluster NGC 6624. iGuver et all (1201 Obi) uses a boxcar 



distri bution from 6.8 kpc t o 9.6 kpc, to reflect two distance measurements of (7.6 ± 0.4) kpc (IKuulkers et al. 



2003h and (8.4 ± 0.6) kpc dValenti et al.ll2007h . We employ a Gaussian distribution centered at 8.2 kpc with 
CT]) = 0.7 kpc since th e 95% c onfidence regions for this Gaus sian is the same as the range suggested by the 
boxcar in lGuver et all (|2010b|). A more recent determination (|Dotter et al.ll2010h gives an apparent distance 



of 10.2 kpc, which when corrected for extinction gives 8.1 kpc, consistent with our value described above. 
Touchdown flux: Five of the bursts have a measured touc hdown flux: FTnoo/ (10~^ ergs cm~^ s~^) 



533±0.27. 5.65 ±0.20. 5.12 + 0.15. 5.24±0.19. and 5.42 ±0.16 (lGuveretal.ll2010bh. The fit in buveretal. 
(l2010br) gives (5.31 ± 0.10) x 10"*^ ergscm"^ s"^, which we use here. 



Hydrogen mass frac tion: This object is likely an ultra-compact binary with a H-poor donor (|King & Watson 



19861 : IStellaet al.lll987h . Altho ugh evolutionary models d o not exclude the possibility that the envelope 
could contain some H (X < 0.1: IPodsiadlowski et all 2002r) . a comparison of burst properties with theoreti- 
cal ignition models suggests H-poor fuel (ICummingll2003b . We fix X = for this source. 
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Table 9. Most probable values for masses and radii for neutron stars constrained to lie on one mass versus 

radius curve. 



Object 


M(Mo) 


/?(km) 


M{Mq) 


7? (km) 




rph 


= R 


fph 


»/? 


4U 1608-522 


1 52+°-22 


1 1 04+0-53 
^^■" -1.50 


1 64+0-34 


11.82_+o.42 


EXO 1745-248 


1 55+0-12 

^"'■'-0.36 


10 91+0-86 

^"■^ -0.65 


, 04+0.450 
^-' -0.28 


11.82:°:^^ 


4U 1820-30 


1 57+"" 

^"^ -0.15 


10 91+0.39 


1 57+0-37 

^""-0.31 


11.82!°:^2 


M13 


1 48+*'-2i 


11.04!}oo 


0.90l!0-28 


12.2l!0.iB 


a»Cen 


1 4O+0.26 


11.18!}i4 


0.994+°-^} 


12 09+*'-27 
^^•"^-0.66 


X7 


o-832:i:i^i 


13 25+1-3' 


1 98+0-10 

^•^^'-0.36 


1 1 O+0.95 
^^•■'-1.03 
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